clear all
%%%%% load the data for validation
load('\\enac1files\change\common\Teaching\Ecohydrological modeling\Exercises\T&C modelling\Exercise 1D\Swiss_Case_Studies\Pre_process\Res_CH-Lae.mat')
k=1; %% GPP type

Date_obs=Date;
%%% set up observed data
LE_obs=LE; %% latent heat
H_obs=H; %% sensible heat
Rn_obs=Rn; %% Net radiation
%%% GPP
GPP_all=[GPP1,GPP2,GPP3,GPP4];
GPP_obs=GPP_all(:,k);

%%% load simulation results
load('C:\Users\tlian\Desktop\Ecohydrology teaching\Exercise 1\Swiss_Case_Studies\Swiss_Case_Studies\results\LAEGEREN.mat','NPP_H','NPP_L','RA_H','RA_L','Rn','H','QE','Ccrown','dateNum1','dateNum2')
%%% daily GPP calculation
GPP_H=(NPP_H+RA_H);
GPP_L=(NPP_L+RA_L);
GPP_sim=(GPP_H+GPP_L)*Ccrown';
LE_sim=QE;


%%% time index
index1 = find(Date_obs >= dateNum1, 1, 'first');
index2 = find(Date_obs >= dateNum2, 1, 'first');
x1=index1;
x2= index2;

LE_obs=LE_obs(x1:x2); %% latent heat
H_obs=H_obs(x1:x2); %% sensible heat
Rn_obs=Rn_obs(x1:x2); %% Net radiation
GPP_obs=GPP_obs(x1:x2);
GPP_obs_daily = mean(reshape(GPP_obs, 24, []))';
%%% unit change from µmol CO₂/m²/s to gC/m²/day
GPP_obs_daily=GPP_obs_daily*12*10^-6*24*3600; % gC/m²/day

Date=Date(x1:x2);

%%% latent heat
figure(1),
plot(Date,LE_obs)
hold on
plot(Date,LE_sim)
ylabel('Latent heat [W/m^2]')
datetick('x', 'dd/mm/yyyy', 'keeplimits', 'keepticks');
set(gca,'FontSize',14)
legend('OBS','SIM')

figure(2)
x = LE_obs(~isnan(LE_obs));
y = LE_sim(~isnan(LE_obs));
axis equal
scatter(x,y)
hold on
ylabel('SIM [W/m^2]')
xlabel('OBS [W/m^2]')
plot([0, 800], [0, 800], 'k--')
set(gca,'FontSize',20)
p = polyfit(x, y, 1);
y_fit = polyval(p, x);
SS_res = sum((y - y_fit).^2);
SS_tot = sum((y - mean(y)).^2);
R2 = 1 - SS_res/SS_tot;
xlim([0 800])
ylim([0 800])
text(0.1, 0.9, ['R^2 = ', num2str(R2, '%.2f')], 'Units', 'normalized', 'FontSize', 20);
box on


time_plot=[Date(1):Date(end)];
figure(3),
plot(time_plot,GPP_obs_daily)
hold on
plot(time_plot,GPP_sim(2:end))
xlim([Date(1) Date(end)])
ylabel('Daily GPP [gC m^{-2} day^{-1}]')
datetick('x', 'dd/mm/yyyy', 'keeplimits', 'keepticks');
set(gca,'FontSize',14)
legend('OBS','SIM')

figure(4)
x0 = GPP_obs_daily;
y0 = GPP_sim(2:end);
x=x0(~isnan(x0)); %% whether exclude negative values. or we keep both
y=y0(~isnan(x0));
axis equal
scatter(x,y)
hold on
ylabel('SIM [gC / m^2 day]')
xlabel('OBS [gC / m^2 day]')
plot([0, 30], [0, 30], 'k--')
set(gca,'FontSize',20)
p = polyfit(x, y, 1);
y_fit = polyval(p, x);
SS_res = sum((y - y_fit).^2);
SS_tot = sum((y - mean(y)).^2);
R2 = 1 - SS_res/SS_tot;
xlim([0 30])
ylim([0 30])
text(0.1, 0.9, ['R^2 = ', num2str(R2, '%.2f')], 'Units', 'normalized', 'FontSize', 20);
box on

%%% sensible heat
figure(5),
plot(Date,H_obs)
hold on
plot(Date,H)
xlim([dateNum1 dateNum2])
ylabel('Sensible heat [W/m^2]')
datetick('x', 'dd/mm/yyyy', 'keeplimits', 'keepticks');
set(gca,'FontSize',14)
legend('OBS','SIM')

figure(6),
x = H_obs(~isnan(H_obs));
y = H(~isnan(H_obs));
axis equal
scatter(x,y)
hold on
ylabel('SIM [W/m^2]')
xlabel('OBS [W/m^2]')
plot([0, 800], [0, 800], 'k--')
set(gca,'FontSize',20)
p = polyfit(x, y, 1);
y_fit = polyval(p, x);
SS_res = sum((y - y_fit).^2);
SS_tot = sum((y - mean(y)).^2);
R2 = 1 - SS_res/SS_tot;
xlim([0 800])
ylim([0 800])
text(0.1, 0.9, ['R^2 = ', num2str(R2, '%.2f')], 'Units', 'normalized', 'FontSize', 20);
box on


%%%% net radiation heat
figure(7),
plot(Date,Rn_obs)
hold on
plot(Date,Rn)
xlim([dateNum1 dateNum2])
ylabel('Net radiation [W/m^2]')
datetick('x', 'dd/mm/yyyy', 'keeplimits', 'keepticks');
set(gca,'FontSize',14)
legend('OBS','SIM')

figure(8)
x = Rn_obs(~isnan(Rn_obs));
y = Rn(~isnan(Rn_obs));
axis equal
scatter(x,y)
hold on
ylabel('SIM [W/m^2]')
xlabel('OBS [W/m^2]')
plot([0, 800], [0, 800], 'k--')
set(gca,'FontSize',20)
p = polyfit(x, y, 1);
y_fit = polyval(p, x);
SS_res = sum((y - y_fit).^2);
SS_tot = sum((y - mean(y)).^2);
R2 = 1 - SS_res/SS_tot;
xlim([0 800])
ylim([0 800])
text(0.1, 0.9, ['R^2 = ', num2str(R2, '%.2f')], 'Units', 'normalized', 'FontSize', 20);
box on
